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Abstract 

^^ I Calculations of observables in quantum chromodynamics are typically per- 

formed using a method that combines numerical integrations over the mo- 

^ ■ menta of final state particles with analytical integrations over the momenta 

of virtual particles. I describe the most important steps of a method for 
performing all of the integrations numerically. 
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I. INTRODUCTION 

This paper concerns a method, which was introduced in [^ , for performing perturbative 
calculations in quantum chromodynamics (QCD) and other quantum field theories. The 
method is intended for calculations of quantities in which one measures something about 
the hadronic final state produced in a collision and in which the observable is infrared 
safe - that is, insensitive to long-distance effects. Examples include jet cross sections in 
hadron-hadron and lepton-hadron scattering and in e~^e~ — > hadrons. There have been 
many calculations of this kind carried out at next-to-leading order in perturbation theory. 
These calculations are based on a method introduced by Ellis, Ross, and Terrano in the 
context of e^e'^ —>■ hadrons. Stated in the simplest terms, the Ellis-Ross- Terrano method 
is to do some integrations over momenta ii analytically, others numerically. In the method 
discussed here, one does all of these integrations numerically. Evidently, if one performs all 
of the integrations numerically, one gains flexibility to quite easily modify the integrand. 
There may be other advantages, as well as some disadvantages, to the numerical integration 
method compared to the numerical/analytical method. 

In this paper, I address only the process e~^e~ -^ hadrons. I discuss three-jet-like infrared 
safe observables at next-to- leading order, that is order a"^. Examples of such observables 
include the thrust distribution and the fraction of events that have three jets. 

The main techniques of the numerical integration method for e~^e~ -^ hadrons were pre- 
sented briefly in [^ . The principle purpose of this paper is to explain in detail some of the 
most important of these techniques. In the numerical/analytical method, one has to work 
hard to implement the cancellation of "collinear" and "soft" divergences that occur in the 
integrations. In the numerical method, as we will see, this cancellation happens automat- 
ically. On the other hand, in the completely numerical method one has the complication 
of having to deform some of the integration contours into the complex plane. We will see 
how to do this deformation. In both the numerical/analytical method and the completely 
numerical method, one must arrange that the density of integration points is singular near 
a soft gluon singularity of the integrand (even after cancellations). However, the precise 
behavior of the densities needed in the two cases is different. We will see what is needed in 
the numerical method. 

These techniques are presented in Sees. |n[-|VI|. They are illustrated in Sec. |V11| with a 
numerical example. Although a full understanding of the example requires the preceding 



sections, the reader may want to look briefly at Sec. |VII| before starting on Sees. |II|-|VI|. A 
brief summary of techniques not presented in detail in this paper is given in Sec. |V111| . 

In Ijll], I presented results from a concrete implementation of the numerical method in 
computer code. Since then, one logical error in the code has been discovered and fixed and 
the performance of the program has been improved. Results from the improved code are 
presented in Sec. ffX. 



Let us begin with a precise statement of the problem. We consider an observable such as 
a particular moment of the thrust distribution. The observable can be expanded in powers 
of aJiT, 

a = ^aW, aWoc(a,/7rr. (1) 



The order a^ contribution has the form 
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Here the daj^^ are the order a^ contributions to the parton level cross section, calculated 
with zero quark masses. Each contains momentum and energy conserving delta functions. 
The dalf^ include ultraviolet renormalization in the MS scheme. The functions S describe 
the measurable quantity to be calculated. We wish to calculate a "three-jet-like" quantity. 
That is, ^2 = 0. The normalization is such that 5„ = 1 for n = 2,3,4 would give the 
order a^ perturbative contribution the the total cross section. There are, of course, infrared 
divergences associated with Eq. (H). For now, we may simply suppose that an infrared cutoff 
has been supplied. 

The measurement, as specified by the functions iS„, is to be infrared safe, as described 
in Ref. |^: the Sn are smooth functions of the parton momenta and 

Sn+i{ki, ..., Xkn, (1 - A)4) = Snih, . . . , 4) (3) 

for < A < 1. That is, coUinear splittings and soft particles do not affect the measurement. 
It is convenient to calculate a quantity that is dimensionless. Let the functions Sn be 
dimensionless and eliminate the remaining dimensionality in the problem by dividing by ctq, 
the total e~^e~ cross section at the Born level. Let us also remove the factor of (as/vr)^. 
Thus, we calculate 

Our problem is thus to calculate X. Let us now see how to set up this problem in a 
convenient form. We note that T is a function of the cm. energy ^/s and the MS renormal- 
ization scale yU. We will choose /i to be proportional to a/s: fi = Ayy^. Then X depends 
on A. But, because it is dimensionless, it is independent of ^/s. This allows us to write 

/■oo 

T= / d^sh{^s)Z{Auy,^s), (5) 

JO 

where h is any function with 

/ d^Ts K^fs) = I. (6) 

Jo 

The quantity X can be expressed in terms of cut Feynman diagrams, as in Fig. ^ The 
dots where the parton lines cross the cut represent the function iS„(A;i, . . . , k^). Each diagram 
is a three loop diagram, so we have integrations over loop momenta £i, ^2 ^^^d £3. We first 




FIG. 1. Two cuts of one of the Feynman diagrams that contribute to e^e -^ hadrons. 

perform the energy integrations. For the graphs in which four parton hnes cross the cut, 
there are four mass-shell delta functions S{k'j). These delta functions eliminate the three 
energy integrals over i^, £2, and £3 as well as the integral (^) over ^/s. For the graphs in 
which three parton lines cross the cut, we can eliminate the integration over -y/s and two of 
the £j integrals. One integral over the energy E in the virtual loop remains. We perform 
this integration by closing the integration contour in the lower half E plane. This gives a 
sum of terms obtained from the original integrand by some algebraic substitutions, as we 
will see in the following sections. Having performed the energy integrations, we are left with 
an integral of the form 

X= f di, di2 di, Y. E 9{G, c- fi, /2, k). (7) 

•^ G C 

Here there is a sum over graphs G (of which one is shown in Fig. |^) and there is a sum over 
the possible cuts of a given graph. 

The problem of calculating 1 is now set up in a convenient form for calculation. If 
we were using the Ellis-Ross- Terrano method, we would calculate some of the integrals in 
Eq. (|^) numerically and some analytically. In the method described here, we first perform 
certain contour deformations, then calculate all of the integrals by Monte Carlo numerical 
integration. In the following sections, we will learn the main techniques for performing the 
integrations in Eq. (|^. We will do this by studying a simple model problem that will enable 
us to see the essential features of the numerical method with as few extraneous difficulties 
as possible. 



II. A SIMPLIFIED MODEL 

In the following sections, we consider a simplified model in which all complications that 
are not needed for a first understanding of the numerical method are stripped away. The 
model is represented by the graph shown in Fig. |^. There are contributions from all of 
the two and three parton cuts of this diagram, as shown in Fig. ^ Since QCD numerator 
functions do not play a major role, we consider this graph in 0^ theory. Thus, also, we can 
avoid the complications of ultraviolet renormalization. We consider the incoming momentum 
q to be fixed and nonzero. We calculate the integral of the graph over the incoming energy 
g°. This is analogous to the technical trick of integrating over ^/s in the full three loop QCD 
calculation (see Sec. |ID and serves to provide three energy integrations to perform against 
three mass-shell delta functions for the three-parton cuts. 

We need a nontrivial measurement function S. As an example, we choose to measure 
the transverse energy in the final state normalized to the total energy: 




FIG. 2. Diagram for a simple calculation. All two and three parton cuts of this diagram in cfi^ 
theory are used, with a measurement function that gives the average transverse energy in the final 
state. 



S2{Kk2) = {\kT,l\ + \kT^2\)/{\kl\ + \k2\) 

53(^1,^2,^3) = {\kT,i\ + \kT,2\ + \kT,-i\)/{\ki\ + 1^21 + \kz\ 



(8) 



where kx^j is the part of the momentum kj of the jth final state particle that is orthogonal 
to q. 

There are two loops in our diagram. We choose the independent loop momenta to be C.2 
and ^4. The other momenta are understood to be expressed in terms of I21 ^4, ^"^^ Q^- 







FIG. 3. The two and three parton cuts of the simple ip'^ diagram. 
Thus the example integral that we seek to calculate is 

^^^ fdf f d% f d% ^^, 



2 J 271 J (27r)4 J (27r) 



(9) 



Here g is the coupling, 1/2 is the statistical factor for this graph, and the integrand W 
consists of four parts, one for each of the cuts in Fig. ^: 



W = Wa + Wfe + We + W, 



(10) 



where 

>V„ = iS,{UM) ^^ ^^ ^j^ (27r)A(£4) (27r)A(4), 

-t]^ ~r t-c -02 ~r t-c ^q ~r "^c 

W, = -iS,{l, -4) (27r)A(£0 -^ (27r)A(-£3) j^^ j^. 

We = sSx. -^2,4) (2vr)A(£i) (27r)A(-4) ^^ (27r)A(4), 

€3 €4 

W, = 53(4, ^2, -4) 4 (27r) A(£2) (27r) A(-4) (27r) A(£4) ^. (11) 

Here we have used the notation 

A(A;) = 5(A;2)^(A;0). (12) 



III. THE INTEGRATION OVER ENERGIES 

We begin by performing the integrals over the energies in Eq. (^. In the case of three 
partons in the final state, the three delta functions eliminate the three integrations. In 
the case of two partons in the final state, the two delta functions eliminate two of the 
energy integrations. This leaves one integration over the energy that circulates around the 
virtual loop. There are three poles in the upper half plane and three in the lower half 
plane. Closing the contour in one half plane or the other gives three contributions. Each of 
these contributions corresponds to putting one of the particles in the loop on shell. Thus 
altogether there are eight contributions to X, as indicated in Fig. |^. We write X as 

where the integrand Q has eight parts: 

Q = Qa\^ Ga2 + Ga3 + Qb5 + Qb2 + Qbi + Gc + Qd- (14) 

The contributions to Q are 

Qal = S2{U," " 
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FIG. 4. The eight contributions to the sample diagram after performing the energy integrations. 
The hne through a propagator in a loop indicates that this propagator is put on shell, with positive 
energy flowing in the direction of the arrow. The direction for positive energy flow around the loop 
depends on whether the contour over loop energy is closed in the upper or the lower half plane. 

So far, the operations that we have performed have been purely algebraic. They are 
evidently of a sort that can be easily implemented in a computer program in an automatic 
fashion. We are left with an integral over the loop momenta £2 and £4. We seek to per- 
form this integration numerically. However, the integrand Q has singularities, so it is not 
completely self-evident how to proceed. It is to this question that we now turn. 

IV. CANCELLATION OF SINGULARITIES 



In this section, we discuss the cancellation of singularities in a numerical calculation of 
the integral in Eq. ([13|) . 

Let us concentrate to begin with on the cut shown in Fig. 0(a). Then there is a virtual 
loop consisting of the propagators with momentum labels ii, £2 and 4- Recall that we are 
taking £2, and 4 as the independent loop momenta. Put the integration over 4 inside the 
integration over 4- Then we can consider 4 as fixed while 4 varies. Fig. |^ illustrates the 
space of the loop momentum £2 for a particular choice of q and at a particular point in the 
integration over 4- The origin of coordinates is at the point labeled 4 = 0. The vector 4 
is indicated as an arrow with its head at £2 = 0. Then the point £1 = is at the tail of this 



vector, as indicated. The vector £5 = g — ^4 is indicated as an arrow with its tail at £2 = 0. 
Then the point £3 = is at the head of this vector, as indicated. Finally, the vector q is 
indicated as an arrow with its tail at i\ = 0. 



1=0 




L=0 



L=0 



FIG. 5. Space of loop momentum 
sentative choice of q, £4, and £5 = q — 



for the virtual loop in the graph of Fig. 0(a) for a repre- 



Where are the singularities of the integrand for our graph? 

There is, first of all, a singularity when the momentum of any propagator vanishes since 
there is always a contribution in which that propagator is put on-shell, with a singularity 
1/(2|£|). Since an integration J di/{2\i\) is convergent in the infrared by two powers, these 
singularities do not cause much difficulty. We simply have to choose a density of points 



with a matching l/\i\ singularity, as described later in Sec. ^ We do not discuss these 
singularities further in this section. 

The singularities of concern to us here are 

1) A coUinear singularity at £2 = —xd^ with < x < 1. 

2) A collinear singularity at £2 = ^(^^ with < a; < 1. 

3) A soft singularity at £2 = 0. 

4) A scattering singularity at |£i| + l^sl = \(ii\ + l^sj. 

The locations of these singularities are indicated in Fig. 1^. 



I^J + l^sl = 1^41 + 14 




-to — XL^ 



t=0 



FIG. 6. Locations of singularities of Qa- 



A. The collinear singularities 

In this subsection, we examine the coUinear singularity at £2 = —^£4 with < x < 1. 
The principles that we discover for this case will hold for the other collinear singularities as 
well. 

The terms Qai and Q^ in the integrand Q, Eq. ([T^), are singular along the line £2 = —xi^, 



< X < 1. In order to examine this singularity, let us write Qai and Qc as given in Eq. (15) 
in the form 



2 'JZ{Ei,Ei'\E,,e2,h) Ss{h,-i2,q-Q. (17) 



21^2 + ^41 mi [Ei-E. 



2 i *-4 



Here the first factors exhibit the denominators for the three propagators that carry collinear 
momenta at the singularity, TZ denotes the rest of the Feynman graph, and the S functions 
are the measurement functions for the final state particles. The functions TZ depend on 
the loop momenta I2 and £4 and on three loop energies, which we take to be Ei = C\, 
E2 = £2 and -E5 = £5. The energies are determined by the on-shell delta functions for the 
two contributions. For Ei and E^, the values are the same for the two contributions: 

^1 = 1/2 + ^41, 

E, = \q~h\. (18) 

For E2, the values are different: 

Et^ = \i2 + i4\-\i4\, 
E^^ = -u. (19) 

In order to examine the behavior of Qai and Qc near the singularity, let 

fa = -xii + /t, (20) 

where £t ■ ^4 = 0. The singularity is at fj- — > 0. 

In Qai the denominator (-E2 ) ~ ^2 vanishes at £t — ^ 0: 

(Er)'-/| = --^(l + 0(/|)). (21) 



Thus there is a 1/^^ singularity which would give a logarithmically divergent result for the 
integral of Qai alone. Altogether, the denominator factors for Qai are 

' ' ' ' i (1 + 0(4)) • (22) 



2|^2+^4| (Ef)) -f| 2\U\ 



t2 02 



Let us now look at the denominator factors for Qc- The denominator (Ei — E. 
takes the form 



W^' ^2 



El - Ei'^^Y - flV = — ^ (l + 0(4)) , (23) 



x{l — x) 
so that the denominator factors together take the form 



214 + 41 2|4| (E,~Ei'T-ii 4£| £| 



^^^ = i i (l + ^(^1)) • (24) 



Again, we have a l/ij^ singularity. 

Note, however, that the denominator factors in Eqs. ( p2D and ( plD are equal except for 
their sign, up to corrections that are not singular as £^ — > 0. Thus if the remaining factors 
TZ and S were exactly the same for Qai and Qc there would be no singularity in their sum. 

We thus need to explore the matching of TZ and S. The two versions of TZ are the same 
functions with the same arguments except for the fact that E2 7^ E2 ■ However, 

Ef) = Et^ + 0{i^). (25) 

Thus 

7^(El, E^^\ E,, 4, h) = n{E^, Ef \ ^5, h, h) + 0{ll). (26) 

For the functions S used in our example, we have 

53((1 - x)U + 4, xU - 4, q- 4) = 52(4, g"- 4) + 0{il). (27) 

Using these matching equations we find that 

QaX^Qc = 0{\) (28) 

as £r ~^ 0. There is no coUinear singularity in Q. 

How general is this result? First of all note that, in the part of the argument not involving 
the measurement functions iS, we used only the explicit structure of the denominators for 
three propagators that meet at a vertex. In the limit in which the momenta carried on these 
propagators become coUinear, there is a cancellation of the coUinear singularity arising from 
these denominators. The three propagators can be part of a much larger graph, and there 
can be non-trivial numerator factors, as in QCD. All of the other factors can be lumped 
into a function IZ and treated as above. Thus this cancellation works in QCD as well as 0^ 
theory and it works for cut graphs with at most one virtual loop at any order of perturbation 
theory. 

As for the measurement functions, in general we need to consider the difference between 
the measurement functions with n and n + 1 particles in the final state, 

F{Jt) = <Sn+i{ki, ■ • ■ , K-i,xkn - It, (1 - x)kn + /t) - Sn{ki, ■ ■ ■ , 4-1, kn). (29) 

10 



Assuming that F is an analytic function of It-, it will have an expansion around ^y = of 
the form 

F(Jt) = a + bi-i^ + Ci/^4 + ■ ■ ■ ■ (30) 

Infrared safety requires that a = 0. If foj 7^ then F vanishes on a surface that intersects 
the point It = 0. Measurement functions S with this property would define an infrared 
safe measurement, but I do not know of any example in common use. More typically, F is 
non-zero in a neighborhood oi ix = while vanishing at £t = 0. Then both a and the bi 
must vanish and the Cy should be a positive definite (or negative definite) matrix. Thus, for 
typical measurement functions, 

F(/t) = 0(4) (31) 

as £t ~^ 0. Then the integrand does not have coUinear singularities. 

For an atypical measurement function with fej 7^ 0, one would be left with an integrable 
singularity of the form h-lx/^T- The current version of the computer code has a mechanism 
to deal with this contingency, but I do not discuss it here since I know of no case in which 
it is needed. 



B. The soft singularities 

— * 

In this subsection, we examine the soft singularity at £2 = 0. 

Let us concentrate to begin with on the cut graph shown in Fig. ^(a). When we perform 
the integration over the energy circulating in the virtual loop, there is a contribution from 
the term in which the propagator carrying momentum ^^ is put on shell, as in Fig. §(al). 
This contribution is Qai in Eq. (0). Let us examine this contribution in the limit £2 -^ 0. 
Expanding in powers of I21 we have 

£? = |/l| = 1/4 + U = \^ + \UU2 ■ M4 + ■ ■ ■ , (32) 

where we adopt the notation 

uj = ij/\Ij\. (33) 

Then 



and 



Thus 

1 



e^ = -\E,\^l-[u2.u,f] + --- (34) 



^3 = 2|4| I/2I U2 ■ (m5 - M4) + ■ ■ ■ . (35) 



Gal ~ S2 



2|4| \kV [1 - {U2 ■ u^f] 2%\ 1^21 «2 ■ (m5 - M4) + It 2141 2|4| 
52 1 1 1 



16 |f4|2 |4|2 fy^ 1 - (m2 ■ U^Y U2 ■ {U5 - U4) + is' 
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(36) 



We proceed in this fashion to evaluate the contribution corresponding to Fig. §(a2). Then 
we evaluate the contribution of Fig. §(a3), but we find that this contribution is not singular 
as ^2 ^ 0. Adding the three contributions, we obtain the net integrand for the cut graph of 
Fig. 1(a) in the soft hmit fa ^ 0: 

" ^ 32 |f4|2 l/sP ifaP l + U2-Uil-U2-u^U2-{u^-Ui) + ie' 

Some comments are in order here. First, we have included the leading term, with a 
l/l^al^ singularity, and dropped less singular terms. If we decompose the integration over 
£2 into ^dVt2 J |^2p'^|^2|, then a 1/|^2|^ singularity produces a logarithmic divergence in the 
integration over 1^2 1- The less singular terms will lead to a finite integration over 1^2 1, 
although the integration JdQ over the angles U2 can still be divergent. There are, in fact, 
singularities in the angular integration. The factor 1/[1 — U2 ■ U5] is singular when £2 is 
collinear with £5, while the factor 1/[1 + ^2 ■ Ua] is singular when —£2 is coUinear with £4. 
These singularities produce logarithmically divergent integrations over U2. However, the 
analysis of the previous subsection shows that the collinear singularities cancel among the 
cuts of our graph. There is also a singularity on the plane U2 ■ (ms — M4) = 0. This is the 
scattering singularity on the ellipse \£i\ + l^al = \£4\ + \£5\. This ellipse passes through the 
point £2 = and the plane tangent to the elhpse at this point is the plane U2 ■ [u^ — u^) = 0. 
I will have more to say about this singularity later. Here we note simply that it comes with 
an ie prescription, which has been preserved in Eq. (|37|). 

We now consider the cut graph shown in Fig. ^(b). Again, there are three contributions 
to consider, corresponding to the diagrams (65), (62) and (64) in Fig. |[ Adding the three 
contributions, we obtain the net integrand for the cut graph of Fig. ^(b) in the soft limit 
£2 -^ 0: 

^ ^ '^2 _J- 1 1 2 + M2 • (m5 - ^4) /ggv 

'' ~ 32 |f4|2 l/sP \£^\3 l-U2-U^l + U2-U5U2-{u5-U^) + ie' 



As in Eq. (pT]), there are a scattering singularity and two collinear singularities. However, 
the signs that indicate the location of the collinear singularities are reversed compared to 
Eq. (pSl). If we add Qa and Qb we obtain 

Q ,Q -^2 1 I + {U2 ■ U4){U2 ■ U5) /oq^ 

'^ ' ^ 16 I/4P I/5P \h\' [1 - (^2 ■ u,ni - iu2 ■ u,y] ■ ^ ^ 

Thus, the overall 1/|^2|'^ singularity remains and the collinear singularities remain, but the 
scattering singularities cancel in the soft limit, £2 — > 0, between the two cuts that leave 
virtual subgraphs. 

There are two more cut graphs to consider. The graph shown in Fig. ^(c) gives 

^' ' ' 1 (40) 



32 |£4p I4P 1^2^ [1 + ^2 ■ ua] [1 + U2- ui 
The graph shown in Fig. |^(d) gives 
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g, S^-^ -1 I I . (41) 

32 1^4^ I4P |^2|3 [1 - M2 ■ M4] [1-U2- M5] 

Adding these together, we find 

'^ ' ^ 16 W I4P I/2P [1 - (^2 ■ «4)2][1 - {U2 ■ U,Y] ■ ^ ' 

We note that when we add the contributions of the cuts which leave virtual subgraphs 
to the contributions of the cuts which have no virtual subgraphs, the leading soft singularity 
cancels: 

Ga + Qb + Gc + Qd-^ 0. (43) 

That is, after cancellation, the overall singularity is at worst proportional to 1/1^2 P- It is 
thus an integrable singularity provided that all of the singularities of the angular integration 
over U2 cause no problems. 

The cancellation of the leading soft singularity is built into the structure of Feynman 
diagrams, so that we do not have to do anything special to make it happen. However, 
there is a certain subtlety in arranging for the singularities in the angular integrations to 
be convergent in a Monte Carlo integration. Thus, we will return to the cancellation of the 
soft singularity after we have discussed contour deformations in the following section. 



V. THE SCATTERING SINGULARITY AND CONTOUR DEFORMATION 

Consider the contribution from Fig. §(al), as given in Eq. ([l5|) . There is a factor 
1 1 



(|/i| - \U - \U? - ?i + ^e (141 + \U + 141 - |4|)(|/4| + 141 - 141 - 141 + ^e) ' 



(44) 



which has a singularity when |£i| + |4| = |4| + l^sl- In an analysis using time-ordered 
perturbation theory, the singular factor emerges from the energy denominator associated 
with the intermediate state consisting of partons 1 and 3, 

Ef - ^(4) + le, (45) 

where Ep = |4| + 14 1 and 

^(4) = 141 + 141 = 14 + 41 + I - 4 + 4|. (46) 

Thus the singularity appears when the momenta are right for particles 1 and 3 to be on-shell 
and scatter to produce the final state particles 4 and 5. 

The contribution from Fig. H(b4) has a scattering singularity at the same place as that 
from the cut diagram (al). However, these singularities do not cancel in general because 
the functions 52(4, -^5) and iS2(£i,4) do not match. We thus have a problem if we would 
like to perform the integration numerically. 
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We notice, however, that the singularity is protected by an ie prescription. The ie in the 
denominator tells us what to do in an analytic calculation and it also tells us what to do in 
a numerical calculation: we need to deform the integration contour. 

We are integrating over a loop momentum £2- Let us replace £2 by a complex momentum 
£2,0 = £2 + m, where k is a function, which remains to be determined, of £2- Then as we 
integrate over the real vector £2 , we are integrating over a contour in the space of the complex 
vector £2,0- When we deform the original contour £2,0 = £2 to the new contour £2,c = £2 + i^, 
the integral does not change provided that we do not cross any points where the integrand 
is singular and provided that we include a jacobian 

^(^-^)-det(^). (47, 

There are some subtleties associated with this; the relevant theorem is proved in the Ap- 
pendix. 

We need to choose /? as a function of £2- Consider first the direction of k. On the 
deformed contour, the energy denominator (|45|) has the form 

Ef - E{£2 + iR) + ie. (48) 

In order to fix the direction of deformation, it is useful to consider what happens when we 
deform the contour just a little way from the real £2 space. For small k, we have 

E{£2 + IK) ^ \£i\ + \U +IK-W, (49) 

where 

^ = 4^ + iL. (50) 

l^il 141 

Thus the energy denominator is Ep — -£^(^2) — ik. ■ w + ie for small K. In order to keep on 
the proper side of the singularity, we want R- w tohe negative. The simplest way to insure 
this is to choose /? in the direction of —w. Thus we choose 

K = -D{£2)w, D{£2)>0. (51) 

Then the singular factor is approximately 



1 



Ef-E{£2) + iD{£2)w' 



(52) 



for a small deformation. For a larger deformation, it not so simple to see that we stay on 
the correct side of the singularity, but it is easy to check numerically. 

The next question is how should we choose D{£2)1 We want D not to be small when £2 is 
near the surface £"(^2) = Ep in order that the integrand not be large there. We want D{£.2) 
not to grow as £| — ^ 00 in order to satisfy the conditions for the theorem that deforming 
the contour does not change value of the integral. Since there is no reason to keep any finite 
contour deformation for large £|, we will simply choose to have £'(£2) — ^ as £| ^ 00. 
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There is another condition that -D(^2) should obey: it should vanish at points where Q 
has collinear and soft singularities. To see why takes some discussion. 

Consider the contributions from three parton cuts, for which there is no virtual loop. 
For these contributions, we do not want to deform the contours. This is because if any 
of the loop momenta were complex then at least one of the momenta of the final state 
particles would be complex. In principle, one could have complex momenta for final state 
particles as long as the measurement functions Sn{ki, . . . , kn) are analytic. However, I have 
in mind applications in which the numerical integration program acts as a subroutine that 
produces "events" with final state particle momenta {ki, . . . , kn} and weights computed by 
the subroutine. Then the events could be the input to, for example, a Monte Carlo program 
that generates parton showers and hadronization. Surely complex momenta for the final 
state particles are not desirable. 

Now recall that there is a cancellation among the contributions Qc from different cuts 
C at points where the Qc have collinear and soft singularities. Evidently, if we deform 
the contour for a contribution with a virtual graph but do not deform the contour for the 
canceling contribution, then the cancellation can be spoiled. We can avoid spoiling the 
cancellation if we make the contours match at the singularity. That is, -D(^2) should vanish 
at the points where the Qc have collinear and soft singularities. 

We also need to determine how fast -D(^2) needs to approach zero as £2 approaches a 
singularity. Since the integration is in a multidimensional complex space, we need an analysis 
that makes use of the multidimensional contour deformation theorem. This analysis is given 
in the Appendix. Here, I present a simpler one dimensional analysis that can serve to clarify 
the issue. 

Consider the following toy integral, 

^"'"^ +/k(.)|. (53) 




-1 + ie 

Here the endpoint singularity at x = plays the role of the collinear or soft singularities. 
The function /y/(x — 1 + ie) plays the role of the integrand for the contribution with a 
virtual subgraph. In this contribution, there is a singularity at a; = 1 that comes with an ie 
prescription. The function fji plays the role of the integrand for the contribution with no 
virtual subgraph. We assume that fv{z) and fniz) are analytic functions. We also assume 
that /v(0) = /r(0), so that the apparent singularity at a; = cancels. 

Now the ie prescription on the singularity at x = 1 tells us that we can deform the 
integration contour into the upper half plane, replacing x hj z = x + iy{x) where y{0) = 
y{xmi,x) = 0. Thus 

' 1 + iy'jx) j fv{x + iy{x)) , • / ^^1 fKA\ 

dx — . . . < -r——T^ + fR[x + iy[x))). (54 

x + iy[x) yx — l + iy[x) J 

Suppose, however, that we want to keep the contour for f^ on the real axis. Then we might 
hope that 1 = 1, where 

~ hm r^dx I '-±^ fvi- + M^)) + /hM 1 (55) 

iiin->o Ja:^j^ [x + iy{x) x — l + iy{x) x J 



X 
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The difference is 

/ - / = lim r^dx |M^ - [1 + .,'(x)]M^i±M^ I . (56) 

x^min^oJx^^^ [ X x + iy[x) ) 

If we note that [fR{z) — /^(0)]/z is an analytic function even at z = and that the integral 
of an analytic function around a closed contour vanishes, we have 



•rnin 



Subtracting these and performing the integral, we have 
i-I = MO) hm r^'^dx j - - 1±M^ I = MO) lini log f 1 + z ^^^^"j . (58) 

We can draw two conclusions. First, as long as y{x) -^ at least as fast as x^ as a: — >■ 0, 
we will realize the cancellation of the x ^ singularity and obtain a finite value for /. 
Second, if we choose y{x) oc x as x ^ 0, I will be finite, but it will not be equal to the 
correct result /. In order to get a result / that is not only finite but also correct, we need 
y{x)/x — i> as X — i> 0. A convenient choice is y{x) oc x^ as x — *> 0. 

We conclude from the multidimensional extension of this analysis, given in the Appendix, 
that as £2 approaches a singularity, -D(^2) should approach zero quadratically with the dis- 
tance to the singularity. 

We now use the qualitative criteria just developed to give a specific choice of deformation. 
We have chosen 

i2,, = e2-tD{e2)w, (59) 

where w, Eq. (0), specifies the direction of deformation. We now specify a deformation 
function .0(^2) that satisfies our criteria. We write D in the form 

D = CG. (60) 

The factor C is designed to insure that the deformation vanishes quadratically near the 
collinear and soft singularities. The factor G is designed to turn the deformation off for 



large £2- These factors are explained below and are defined precisely in Eqs. ( plD and (pTf) 
below. 

First, we discuss the factor C. We want the deformation to vanish at the line £2 = —xi^ 
with < X < 1, where the amplitude has a collinear singularity. (Since £4 = d-i — £2-, this 
line is also £1 = — A£2 with < A < 00.) Define 



K2K1 + K1K2 1^21^1 + |^l|^2 



dl2 = ^ ^ = ^ . (61) 

14 -^2| 141 

This function is zero on the line £2 = — a^4 with < x < 1, and furthermore, it vanishes 
linearly as £2 approaches this line. Similarly, we want the deformation to vanish on the 
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line I2 = 2;^5 with < x < 1, where the amphtude has its other colhnear singularity. The 
function ^23, where 



1^314+1^214 1414 + 1414 



^23 = ^ ^ = ^ , (62) 

14-41 141 

— * —*—#—* 

vanishes linearly as £2 approaches this line. (To see this, use 4 = ^2 — 4-) Let 

d = min(cii2,(i23)- (63) 

Then d vanishes linearly with the distance to either of the collinear singularities. It also 
vanishes linearly with the distance to the soft singularity at £2 = 0. 

Now, we have seen that the deformation should vanish quadratically with the distance 
to any of the singularities. We can achieve this by letting 

where a and (5 are adjustable dimensionless parameters. Note that, for large c?, C[d?) 
approaches a constant. 

Next, we discuss the factor G. We want to ensure that the contour deformation vanishes 
for large ^2- Let us define 

a=|4| + |4|-|g1 (65) 

and 

^=|4| + |4|-|g1. (66) 

Then the singularity that we are avoiding by means of contour deformation is at a = A. We 
can turn the deformation off for a^ Ahy setting 

where 7 is an adjustable dimensionless parameter. 

There is a subsidiary reason for this choice. At the singularity, G = !/[(! + 7)^]- The 
factor 1/A serves to enhance the deformation in the case that 4 and £5 are nearly collinear, 
in which case d is small on the ellipse a = A and the deformation would otherwise be too 
small. 

The reader will note that, while there is a certain uniqueness in defining the direction of 
the deformation in Eq. (^) to be given by the vector w, Eq. (|50|), the normalization D = C G 
with G and G given in Eqs. ( |5^ and ( p7D is rather ad hoc. Within the requirements that the 
deformation should vanish quadratically at the collinear and soft singularities and should 
vanish for large 4, many other choices would be possible. The choice given here is used in 
the current version of the code H]. Surely there is some other choice that is better. 
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VI. THE MONTE CARLO INTEGRATION 

After the contour deformations, we have an integral of the form 

1 = ldiJ2JiC;i) giC;i + tK{C;i)), (68) 



c 



where we use i for the loop momenta collectively, i = {£2, h}- The index C labels the cut, 
a, b, c, or d in Fig. ^ There is a contour deformation that depends on the cut, as specified 
by k(C;£), and there is a corresponding jacobian J'{C;£), Eq. (|47|). Define 



f{i) = ^r£JiC■J)9{C■,i + ^K{C■,i))\. (69) 

We know that I is real, so 

I=ldifii). (70) 

To perform the integration, we use the Monte Carlo method. We choose points £ with a 
density p(i'), with 



dip{i) = 1. (71) 

After choosing N points £1, . . . , ^at, we have an estimate for the integral: 

This is an approximation for the integral in the sense that if we repeat the procedure a lot 
of times the expectation value for Xjy is 

(In) = X. (73) 

The expected r.m.s. error is S, where 

.-((I„-X)^>4/d.^4. (74) 

One can rewrite this as 

where 

i=fdi\f{e)\. (76) 



We see, first of all, that the expected error decreases proportionally to l/yiV. Second, we 
see that the ideal choice of p{i) would be p{i) = \f{i)\/I. 



Of course, it is not possible to choose p in this way. But we know that |/| has singularities 
at places where propagator momenta vanish and we know the structure of these singularities. 
We are not really able to choose p so that |/(^)|/p(^) is a constant, but at least we can choose 
it so that |/(^)|/p(^) is not singular at the singularities of |/(^)|. 

Note that it is easy to combine methods for choosing Monte Carlo points. Suppose that 
we have a recipe for choosing points with a density pi that is singular when one propagator 
momentum vanishes, a recipe for choosing points with a density p2 that is singular when 
another propagator momentum vanishes, and in general recipes for choosing points with 
densities pi with several goals in mind. Then we can devote a fraction Aj of the points to 
the choice with density pi and obtain a net density 

pW = E^^aW- (77) 



A. The density near where a propagator momentum vanishes 

— * 

Let £j be the momentum of one of the propagators in our graph. We have seen that 
when particle J appears in the final state, there is a factor l/|^j| in the integrand. When 
propagator J is part of a virtual loop, the contribution corresponding to this propagator 
being put on shell also contains a factor l/|-^j|. Thus there is a singularity l/|^j| for every 
propagator in the graph. 

The analysis given in the introduction to this section indicates that for each propagator 
J one of the terms pi in the density function should have a singularity that is at least as 
strong as 

P^{^) cx I/I/jI (78) 

as dj — * 0. It is, of course, easy to choose points with a density proportional to l/\aj\-^ as 
£j ^ as long as A < 3. (The limitation on A arises because for A > 3 we would have 
J dij p{i) = oo.) Thus it is easy to arrange that the density of points has the requisite 
singularities. Specifically, we can choose £j with the density 



P(^. 



_2 1 Kq 

2ttK^ 



1 + {\lj\/Ko) 



\l 



(79) 



J\ 



where Kq is a momentum scale determined by the other, previously chosen, loop momenta. 

The singularity when ij ^ can be more severe than l/|-^j|, depending on the structure 

of the graph. Consider first the cases J = 1,3,4,5. Here, the singularities for particular 

cuts, as given in Eq. (IT^), are l/|£jp. However, there is a cancellation after one sums over 



cuts (as for the singularity for £2 -^ 0), leaving a singularity l/|-^j|. 

For J = 2 there is a severe singularity of the form l/|£j|^ for particular contributions to 
Eq. (^). A 1/1^2!^ singularity would not be integrable, but, as we have seen in detail, there 
is a cancellation among the contributions so that only a l/|-^2p singularity is left. However, 
it will not do to simply chose pi{i) oc l/|^2p because there is also a singularity in the space 
of the angles of £2- It is to this subject that we now turn. 
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B. The soft parton singularity 

When two partons can scatter by exchanging a parton before they enter the final state, 
there is a severe singularity as the momentum of the exchanged parton goes to zero. For 
our graph, this happens for £2 ~* 0. In this subsection, we consider the behavior of the 
integrand for small £2 as a function of its magnitude £2 and of its direction U2 = £2/ 1^21 ■ 

The singularity for individual cuts, as given in Eq. (|T5|) , is of the form 1/1^2^ when we 
let 1^21 "*■ with U2 held fixed. This singularity is not integrable. However, as we have seen, 
the leading term cancels when we sum over cuts, leaving a l/|-^2p singularity for 1^21 -^ 
with U2 fixed. 

Let us now recall from Eq. (|^) that, before we deform the integration contour, the 
contribution for small £2 from the cut a of Fig. ^ has, in addition to a factor l/|£2p, a factor 
l/[u2 ■ (^5 ~ W4) + i^]- That is, there is a singularity on a surface in the space of £2 whose 
tangent plane is the plane perpendicular to U5 — U4. We have avoided this singularity by 
deforming the integration contour. However, the deformation vanishes as £2 — > 0. Thus we 
must face the question of what happens to the cancellation near the soft parton singularity 
when the contour deformation is taken into account. 

First, let us recall from Eq. (p9D that for cut a in Fig. ^ the deformation has the form 



£2,c = £2-^D{£2)w, (80) 

where w = ui + u^. For £2 — > 0, 

W r^ Ui — Ms, (81) 

while D has the form 

D{£2)^£l D{u2). (82) 

Here D vanishes for U2 = —u^ and for U2 = u^ and is positive elsewhere. Thus 

£2,0 = I/2I {U2 + 1\£2\D{U2) {u, - u,) + 0(/|)) . (83) 



Substituting £2,c as given above for £2 in Eq. (|37|), we obtain an expression for the 
contribution from cut a to the integrand on the deformed contour near the soft singularity: 

-'^2 1 1 1 2-U2- {U5-U4) 



32 |£4p I4P |4P ^+U2-U^I-U2-U5U2-{U5- U4) + i\£2\D{u2) (m5 - ^4)^ 

There are two cases to consider. First, when 1^21 — * with ^2 fixed, we can drop the second 
term in the last denominator. Then Qa ~ h{u2) / \£2\'^ , where the function h{u2) is the same 
as on the undeformed contour. As we have seen, the leading 1/1^21^ terms cancel when one 
sums over cuts. Thus, as noted earlier, the net integrand behaves like 

G ~ htot{u2)/\£2\' (85) 

when 1^2 1 ^0 with U2 fixed. 
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The second case is more interesting. Consider 1^21 -^0 and U2 ■ [u^ — U4) —>■ with 
U2 ■ {u5 — U4)/|^2| fixed. Then Qa is more singular, Qa oc 1/|£2|^- To see what happens in this 
region, we analyze the contribution from cut b in Fig. ^ in the same fashion. The contour 
deformation for cut b is different from that for cut a, but the deformations match at leading 
order as 1^21 —^ 0. (This is an important feature of the choice of contour deformations.) 



Thus we can use Eq. (|83| ) in Eq. (|38D to obtain 



1 1 1 2 + U2- {U5-U4) 



32 m^ 14 |2 14 P 1 - ^2 • M4 1 + M2 • W5 «2 • («5 - ^4) + 2|4|^(m2) (^5 - U^)^ 

We see that Qb is also proportional to 1/|4|^ in the problematic region. However, since 
M2 ■ M4 ~ M2 ■ M5 in this region, the leading 1/|^2|^ behavior cancels when we add Qb to Qa- 
We are left with the next term, proportional to l/l^l^- 

For the two remaining cuts there is no contour deformation. The contributions from these 
cuts are each proportional to l/|£2p- Calculation shows that there is no further cancellation. 
Thus the net behavior of the integrand is 

Q oc 1/141' (87) 

— * — * 

when |4| — ^ and U2 ■ (ms — U4) — > with U2 ■ {u^ — 'U4)/|4| fixed. 

C. Density near a soft parton singularity 

According to the analysis at the beginning of this section, we should choose a density of 
integration points that has a singularity that is at least as strong as that of |^| near the soft 
singularity at £2 ^ 0. Thus we should choose one of the pi so that 

Pi{t) ex ^— , 141 — > 0, U2 fixed. 



1^2 



P 



Pm oc -J—. \U - 0, ^^^-^f; "^'^ fixed, (88) 

I 2| l'^2| 

with p > 2. 

Specifically, having chosen 4 we can choose the remaining loop momentum £2 with the 
density 



1 1 I KA"^ 1 

P(4 



2txKI 



1+141/^0 



(3_p)-| (5-p)/{3-p) 




vJcos\e) + q/Ki 



(89) 



Here Kq is a momentum scale determined by 4, is the angle between £2 and {u^ — U4), and 

sinh(r) = iro/|4|. (90) 

— * 

It is easy to choose points with this density by first choosing \i2\, then choosing cos{9), and 
finally choosing the corresponding azimuthal angle (p with a uniform density. Accounting 
for the fact that F ex log(|£2|) for £2 -^ 0, we see that p will have a singularity stronger than 
that of Q provided that p > 2. We will see how this works in a numerical example in the 
next section. 
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VII. NUMERICAL EXAMPLE 

In this section, I illustrate the principles developed above by means of a particular 
example. We consider the integral in Eq. (0). We hold £4 fixed and consider the integrand 
as a function of £2- In order to simplify the labelling, I define 

^2 = I (91) 

There is a contribution for each cut C, with C = a,b, c, or d. For each contribution from a 
cut C in which there is a virtual loop, we want to deform the integration contour as discussed 
in Sec. 0. Thus i gets replaced by a complex vector ic = ^ + H^c and we need to supply a 
jacobian Jc{(^)i Eq. (0). Then the integration over I has the form 

fdiY.Jc{i)gcii). (92) 

•^ c 

The functions Qc are the analytic continuations to the deformed contours of the functions 
given in Eq. (|T^). As discussed in Sec. |VJ the quantity that is relevant for the convergence 
of the Monte Carlo integration is the integrand divided by the density of points chosen for 
the integration. In this section, I consider only the integration over i, so I discuss a choice 
for the density of integration points p{i) at a fixed £4 and display plots of the functions 

Fc{i) ^ 4\ ^^^^") ^^(^") (93) 

P(^) 

and F{i) = J2c ^c{^), as well as plots of the deformation and the density. 
For the numerical examples, I choose 



and then take £4 at the point 



Since i^ = a — C-a we have 



g= (3, -0.5,0) (94) 



(2,-1,0). (95) 



:i, 0.5,0). (96) 



The singularities of the functions Qci.^) lie in the plane of £4 and £5, that is the £2 = plane. 
In the plots, I choose £2 = 0, so that we see the effect of the singularities. I plot |/?a|, |k{,|, 
p. Fa, Fh, Fc + Fd and F as functions of d^ and Cy in the domain —2.5 < d^ < 1-0 and 



-1.0 <iy< 2.0. 



Consider first the contour deformation for cut a, ^ ^ ^^ = ^ + 'if^a- I take /?„ = —D w as 
given in Eqs. ( p9l - |6^ ) with a = P = •y = 1. In Fig. |^, I show a graph of \Ka\ versus ix and 



iy. We see that the deformation is not small. I also display in the figure the lines i = —xi4 
with < X < 1 and i = xi^ with < x < 1, where the coUinear singularities for cut a are 
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FIG. 7. Contour deformation for cut a. The absolute value of the imaginary part Kq of ic for 
cut a is plotted against ix and iy at i^ = 0. The deformation Ka vanishes at the two collinear 
singularities for this cut, which are indicated by lines superimposed on the graph. 



\Kh 




FIG. 8. Contour deformation for cut b. The absolute value of the imaginary part Kb of ic for 
cut b is plotted against i^ and £y at iz = 0. The deformation k^ vanishes at the two collinear 
singularities for this cut, which are indicated by lines superimposed on the graph. 



23 



located. We see that, as desired, the deformation vanishes quadratically as d approaches 
these hnes. 

There is a different contour deformation for cut h. The same formulas apply as for cut a 
with the replacements £4 ^h^ £i, £5 <-> —£3, £<->■—£ and with the sign of /? reversed. I show a 
graph of I Kb I versus l^ and ly in Fig. |[ (This figure does not look like Fig. |^ because £ varies 
with £4 held fixed, not with di held fixed as would be needed if we applied the replacement 
^4 <-^ ^1 to Fig. ^) I also display in the figure the lines d = Xi^ with < A and i = —Xi^ 
with < A, where the coUinear singularities for cut b are located. The deformation vanishes 
quadratically as i approaches these lines. 

The jacobian functions J'a{i) and J'b{i) associated with the contour deformations are 
quite unremarkable, so I omit showing them. 




FIG. 9. Density of integration points. The density has three singularities. Only values p < 1 
are shown. 



(97) 



0; 



Consider next the density of integration points. I choose 

p(f) = 0.2pi(fi) + 0.6 p2(/) + 0.2p3(/3), 

as shown in Fig. p. The function pi has a mild singularity as ii —>■ and is given by Eq. 
with ii = £4 — £ and with Kq set equal to 2. The function ps has a mild singularity as £3 
I use the same functional form with £3 = £ — £5 . For p2 , I use the function given in Eq. 
with Kq = 2 and with the power p taken as p = 2.2. Then p2 has a strong l/[|£p'^ log 
singularity as we approach the i = 0. Furthermore, the density of points is largest near the 
plane iy = 0, the plane that is tangent at ^ = to the ellipsoidal surface that (if we turn off 
the deformation) contains the scattering singularity. In order to display the dependence of 
p on angle near i = 0, 1 plot in Fig. |T0| the angle dependent factor in p2, namely the factor 



1^1/^0 



cos-" (6) + iyKi 



(98) 



in Eq. (^), in a region near i = 0. Here cos(6') = iy/\i\. We see that the density of 
integration points is heavily concentrated very near the plane iy = when |^| is small. 
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FIG. 10. Angle dependent factor for the density p2 of integration points associated with the 
soft singularity. For small |£|, points are concentrated near = 0. 

We are now ready to look at the contribution Fa = Ja Gal P to F from cut a. This 
function is displayed in Fig. |TT| with a small rectangle near £ = removed from the graph. 
We see the two coUinear singularities, at £ = —xi^ and at ^ = xl^ with < x < 1. As £ 
approaches one of these singularities, Fa approaches — cxd. 




FIG. 11. Contribution to F from cut a. The domain —1 < i^ < 1, —0.3 < £y < 0.3, which 
contains the soft singularity, has been removed from the plot in order to make the collinear sin- 
gularities visible. The function Fa is singular along lines from (—2,1) to (0,0) and from (0,0) to 
(1,0.5). Only values Fa > —100 are shown. 

— * 

In the standard method for calculating T, we would perform the integration over i 
analytically for the contribution from cut a. Because of the singularities, the integration 
is divergent. However, we can get a finite answer if we regulate the integral by working in 
3 — 2e spatial dimensions. Then the result contains terms proportional to 1/e^ and 1/e as 
well as a remainder that is finite as e ^ 0. 
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What about the contribution to F from cut b, the other cut for which there is a virtual 

"= 

with < A < cx). As with Fa, as i approaches one of these singularities, Fh approaches — cx). 



subgraph? This function is displayed in Fig. |1^ with the same small rectangle near 
removed from the graph. We see the two coUinear singularities, at £ = Xl^ and at £ = 



F. 




-100^ 



FIG. 12. Contribution to F from cut b. The same rectangular domain as in Fig. [T^ has been 
removed from the plot. The function Ff, is singular along lines that extend from (0, 0) to infinity 
in the directions (—2, —1) and (2, —1). Only values -F^ > —100 are shown. 



There are two cuts, c and d, for which there are no virtual subgraphs. In Fig. ^ I show 
the contribution F^ + F^ from these cuts. We see that Fc + F^ approaches +oo at just the 
singularities where Fa and Fh approach — cxd. 

In the standard method for QCD calculations, we would perform the integration over 
i partially numerically for the contribution from cuts c and d. Of course, we would have 
to do something about the collinear and soft singularities, since otherwise we would obtain 
an infinite result. For instance, if we were to use the phase-space slicing method, we would 
slice away a small part of the integration domain near the singularities and calculate its 
contribution analytically in 3 — 2e spatial dimensions in the limit that the region sliced away 
is small. Then we would be left with a numerical integration of Qc + Qd in the remaining 
region (in exactly 3 spatial dimensions). Evidently, the density of points used in the present 
method would not do for this purpose; we would need to expend more points on the region 
near the collinear singularities. 

We see that the standard method for performing the integrations, in which some parts 
of the integrations are performed analytically and some are performed numerically, is, of 
necessity, rather complicated. In the numerical method, we simply combine Qa, Gb, Qc, 
and Qd and integrate numerically. The argument in the preceding sections showed that the 
contributions from the various cuts cancel as one approaches the collinear singularities. This 
is illustrated in Fig. |14|, where I plot Fa + Fh + Fc + F^ versus i^ and iy. We see, first of 
all, that the singular behaviors at the collinear singularities cancel, just as the calculation of 
Sec. 1^ showed. There is also a cancellation at the soft singularity at £ = 0. There is still a 
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FIG. 13. Contribution to F from cuts c and d. The same rectangular domain as in Fig. 11 
has been removed from the plot. The function Fc + F^ is singular along lines that extend from 
(1,0.5) to infinity in the direction (—2,-1) and from (—2,1) to infinity in the direction (2,-1). 
Only values Fc + F^ < 100 are shown. 



singularity in the integrand at £ = 0, but it is integrable and is removed from F by choosing 
a suitable density of points p. Thus F remains less than about 20 everywhere. 

We can see the remnants of the scattering singularity, which is located on an ellipsoidal 
surface that intersects the plane iz = 0. If it were not for the contour deformation, F would 
be singular on this surface, approaching +oo as one approached the surface from one side 
and approaching — oo as one approached the surface from the other side. Since the deformed 
contour avoids the singularity, the singular behavior is removed and we are left with a ridge 



and valley near the ellipsoid. This structure is illustrated in Fig. 15, in which we see a slice 
through the ridge and valley at ir^ = —0.3. Since the amount of deformation vanishes as 
one approaches i = 0, the width of the ridge and valley structure becomes more and more 
narrow as \i\ -^ 0. Recall that the density p of integration points is designed to match this 
increasing narrowness as \i\ — > 0, so that the integration points are concentrated where the 
structure is. 



VIII. OTHER ISSUES 

In the preceding sections, we have seen the most important features of the method of 
numerical integration for one loop QCD calculations. There are other important issues that 
are outside of the scope of this paper. I mention these briefly here. 

First, full QCD has a much more complicated structure than (p^ theory, which was the 
example for this paper. However, the complications of full QCD are in the numerators of 
the expressions representing Feynman diagrams, while the cancellations and the analytic 
structure related to the contour deformation have to do with the denominator structure. 
Thus one can simply generate the numerator structure with computer algebra and carry it 
along. 



27 



F 




FIG. 14. The net function F with the contributions from cuts a, 6, c and d combined. There 
are no coUinear singularities. What remains after cancehation from the soft singularity is removed 
by the density of points in the denominator of F. The remnants of the scattering singularity are 
visible, but there is no actual singularity because of the contour deformation. 



F 




FIG. 15. Structure of F along a slice through Fig. ^ at Ix 
valley structure is very narrow. 



0.03 



-0.3. Note that the ridge and 
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Second, the denominator structure in the example used in this paper is not the only 
denominator structure that one needs to treat. In the QCD calculation for three-jet-like 
quantities in electron-positron annihilation at order a^, there are five possibilities for how a 



virtual subgraph can occur inside an amplitude. The possibilities are indicated in Fig. [16. 
For each possibility there is an entering line representing the virtual photon or Z boson, 
which we take to have zero three momentum, and there are three on-shell lines entering 
the final state. There are graphs of two types, (a) and (b), containing two point virtual 
subgraphs. There are graphs of two types, (c) and (d), containing three point virtual 
subgraphs. There is one type, (e), of graph with a four point virtual subgraph. In structure 
(c), a line with non-zero three momentum enters the three point virtual subgraph and two 
on-shell lines leave. This is the case that we analyzed in the example of this paper. The 
structure of the graph led to the singularity structure depicted in Fig. ^. Amplitudes of types 
(d) and (e) have different singularity structures from that studied here. Case (d) is simpler 
than the case we have studied, while case (e) is somewhat more complicated. However, the 
essential features are those that we have already studied. 




FIG. 16. Ways to insert virtual subdiagrams in Feynman diagrams for e^e -^ 3 partons. 

This leaves virtual self-energy subgraphs. In case (a), there is a self energy subgraph 
on a propagator that enters the final state. This case requires a treatment different from 
that discussed in the previous sections. This is evident because there is a nominal l//c^ 
where fc^ = 0. The treatment required is to represent the virtual self-energy via a dispersion 
relation |T|. In this representation, the subgraph is expressed as an integral over the three- 
momentum in the virtual loop with an integrand that is closely related to the integrand for 
the corresponding cut self-energy graph. The point-by-point cancellation between real and 
virtual graphs is then manifest. It is convenient also to use the dispersive representation for 
the much easier case (b). 

Third, we have to do something about ultraviolet divergences in virtual subgraphs. These 
are easily removed 0] by subtracting an integrand that, in the region of large loop momenta, 
matches the integrand of the divergent subdiagram. The integrand of the subtraction term 
should depend on a mass parameter fj, that serves to make the subtraction term well behaved 
in the infrared. Then, with the aid of a small analytical calculation for each of the one loop 
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divergent subdiagrains that occur in QCD, one can arrange the definition so that this ad 
hoc subtraction has exactly the same effect as MS subtraction with scale parameter /i. 

Fourth, I use Feynman gauge for the gluon field, but then self-energy corrections on 
a gluon propagator require special attention [|l|. The one loop gluon self-energy subgraph, 
71^'^ (k), contains a term proportional to k'^k'^ that contributes quadratic infrared divergences 
0, while the cancellation mechanism that we have studied in this paper takes care of 
logarithmic divergences. This problem can be solved by replacing vr^'^ by P^vr^^P^, where 
P^ = Qa ~ k^ka/k"^, with k = {0,k). The terms added to vr'^'^ are proportional to either 
k^ or fc^ and thus vanish when one sums over different ways of inserting the dressed gluon 
propagator into the remaining subgraph. Since P/^k" = 0, the problematic k'^k" term 
is eliminated. Effectively, this is a change of gauge for dressed gluon propagators from 
Feynman gauge to Coulomb gauge. 

Most of these issues are discussed briefly in ||l[. A quite detailed, but not very peda- 
gogical, treatment can be found in 0. Further analysis of these issues is left for a future 
paper. 

I have also not given a complete presentation of algorithms for choosing integration 
points. As discussed in Sec. |V1C| , the crucial issue is to have the right singularities in the 
density of points near a soft parton singularity of the Feynman diagram. This is not the 
only issue that needs to be addressed in a complete algorithm. Of course, the demonstration 
program [Q has a complete algorithm. However, this algorithm is quite a hodge-podge of 
methods and it seems that a detailed exposition should be reserved for a better and more 
systematic method, which remains to be developed. 

IX. RESULTS AND CONCLUSIONS 

In the preceding sections, we have seen some of the techniques needed for the numerical 
integration method for QCD calculations. Of course, since not all of the techniques have 
been explained, the explanation does not constitute a very convincing argument that such a 
calculation is feasible. A truly exhaustive explanation would help, but an actual computer 
program that demonstrates the techniques is better. Results from such a program were 
presented in [|[]. Since then, I have found and corrected one bug that resulted in errors a 
little bit bigger than 1% and have made some other improvements in the code. The resulting 
code and documentation are available at 0. 

The program is a parton-level event generator. The user is to supply a subroutine that 
calculates how an event with three or four partons in the flnal state contributes to the 
observable to be calculated. The program supplies events, each consisting of a set of parton 
momenta {/ci, /c2, ^3} or {^i, ^2? ^3? ^4}? together with weights w for the events. Then the 
user routine calculates X according to 

1 ^ 
X^-^^,5(fc,). (99) 

The weights used are the real parts of complex weights; the imaginary parts can be dropped 
since we always know in advance that I is real. Thus the weights are both positive and 
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negative. It would, of course, be more convenient to have only positive weights, but one 
can hardly have quantum interference without having negative numbers along with positive 
numbers. 

The first general purpose program for QCD calculation of three-jet-like quantities in 
e^e^ -^ hadrons at order a^ was that of Kunszt and Nason |^. This program uses the 
numerical/analytical method of Ellis, Ross, and Terrano. In Table |, I compare the results 
of the Kunszt-Nason program to those obtained with the numerical method for the a1 
contributions to moments of the thrust distribution, 

^n = , , ,, dt{l- tr^. (100) 

(Jo (as/TTj^ JO at 

In Table Q, I compare the results of the two methods for moments of the y^ut distribution 
for the three jet cross section. To define these quantities, let /sd/cut) be the cross section to 
produce three jets according to the Durham algorithm [§] with resolution parameter ?/cut- 
Let gsiijcut) be the negative of its derivative, 

93{ycu.) = -^f^. (101) 

Then we calculate moments of the a"^ contribution to this quantity, 

In = , \ ,2 f'dycut ivcntrgP. (102) 

(To [as/Try Jo 

In each table, the results for the numerical method are shown with their statistical and 
systematic errors. (The systematic error is estimated by changing the cutoffs that remove 
small regions near the singularities where roundoff errors start to become a problem.) The 
corresponding results for the numerical/analytical method are shown with the statistical er- 
rors as reported by the program. Inspection of the tables shows that there is good agreement 
between the two methods. 

We have explored some of the most important techniques necessary for a QCD calculation 
for three-jet-like quantities in electron-positron annihilation at order a'^ using numerical 
integration throughout the calculation. For the techniques covered, this explanation expands 
on the brief presentation in |]I|. We have also seen that the method works. The older and 
very successful numerical/analytical method for QCD calculations has its complications. The 
numerical method has its own complications, but they are different complications. Thus one 
may expect that the classes of problems for which each of the methods is well adapted may 
be different. There may be some classes of problems for which the natural flexibility of the 
numerical method makes it the more useful method. It remains for the future to explore the 
possibilities. 
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TABLE I. Comparison of results for moments of the thrust distribution, Eq. ( 100 ). The "nu- 
merical" results are from the program 0]. The first error is statistical, the second systematic. The 
"numerical/analytical" results are from the program of Kunszt and Nason [Q] and are given with 
their reported statistical errors. 

n numerical numerical/analytical 

L5 4.127 ±0.008 ±0.025 4.132 ±0.003 

2.0 1.565 ± 0.002 ± 0.007 1.565 ±0.001 

2.5 (6.439 ± 0.010 ± 0.022) x 10^^ (6.440 ± 0.003 ) x 10^^ 

3.0 (2.822 ± 0.005 ± 0.009) x 10^^ (2.822 ± 0.001 ) x 10"^ 

3.5 (1.296 ± 0.002 ± 0.004) x 10"^ (1.296 ± 0.0005) x 10"^ 

4.0 (6.159 ± 0.011 ± 0.016) x 10"^ (6.161 ± 0.002 ) x 10"^ 

4.5 (3.009 ± 0.006 ± 0.007) x 10"^ (3.0IO ± 0.0006) x 10"^ 

5.0 (1.501 ± 0.003 ± 0.003) x 10"^ (1.502 ± 0.0002) x 10"^ 



TABLE II. Comparison of results for moments of the ycut distribution, Eq. ( |102D . The "nu- 
merical" results are from the program [^ . The first error is statistical, the second systematic. The 
"numerical/analytical" results are from the program of Kunszt and Nason [^ and are given with 
their reported statistical errors. 

n numerical numerical/analytical 

L5 (8.442 ± 0.034 ± 0.059) x 10"^ (8.397 ± 0.002 ) x 10"^ 

2.0 (3.106 ± 0.012 ± 0.015) x 10"^ (3.090 ± 0.0004) x 10"^ 

2.5 (1.205 ± 0.005 ± 0.005) x 10"^ (1.200 ± 0.0002) x 10"^ 

3.0 (4.945 ± 0.025 ± 0.019) x 10"^ (4 927 ± q.OOI ) x 10"^ 

3.5 (2.122 ± 0.012 ± 0.008) x 10"^ (2.116 ± 0.0007) x 10^^ 

4.0 (9.430 ± 0.064 ± 0.032) x lO^^ (9.412 ± 0.004 ) x lO'^ 

4.5 (4.304 ± 0.034 ± 0.014) x lO^^ (4.301 ± 0.002 ) x 10"^ 

5.0 (2.008 ± 0.018 ± 0.006) x lO^^ (2.008 ± 0.001 ) x lO^^ 
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APPENDIX: CONTOUR DEFORMATION IN MANY DIMENSIONS 

The calculational method described in this paper makes use of Cauchy's theorem in 
a multi-dimensional complex space. Since this theorem is not proved in most textbooks 
on complex analysis, I provide a proof here, including the special case, needed for our 
application, in which there is a singularity on the integration contour. 

Let f{z) be a function of A^ complex variables z'^ = x'^ + iy^, with /i = 1, . . . , A^, where 
x^ and y^ are real variables. Consider a family of integration contours C{t) labeled by a 
parameter t with < t < 1 and specified by 

z^'{x\---,x^;t) = x^' + iy^'{x\---,x^-t), /x = l,...,iV. (Al) 

Let X{t) be the integral of / over the contour C(t), 

X(t) = j^^^dz f{z) = jdx Aet(^^^\ f{z{x-t)). (A2) 

Suppose that f{z) is analytic in a region that contains the contours. Then we have Cauchy's 
theorem: 

J(1)=X(0). (A3) 

To prove this theorem, we simply prove that dX{t)/dt = 0. Define 

Let Bj^ be the inverse matrix to yl(^. Then 

1 dz"^ dz"'^ 



A2--N 



where e^^'"^^ is the completely antisymmetric tensor with N indices, normalized to e 

(A6) 



1, and e^^...^jY is the same tensor. This has the immediate consequence that 

d 



Also, 



,,.(i^.^detA)^0. 




1 ,m-.., ^^" . 





so 



det A = — e^'-^^e,,...,^ —— ■ ■ ■ -— , (AT) 
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— detA 
dt 

We need one more result: 



d ^ . dK R..„. . . dy^ 



dt 



Bt: det A = i 



dxt'dt " 



Bt'detA. 



so 



df dz"" df 
^x^^ ^x^^ dz" 



df r,, df 



B^ 



Then, using the results 



( |A8|) , and (|A10| ) and an integration by parts, we find 



(A8) 



(A9) 



(AlO) 



d_ 



lit) = fdx — [detA f] 



d X det All 



d X det Ali 



dy"" „„ . , df .dy- 



^x^^^t 



BU + 



dy' 



dz" dt j 
df dy-' 



^x^^^t " ^ " ^x^^ dt 



dy'' ^ d 



i / dx 



/l^WdetA} 



dt dx^^ 



0. 



(All) 



This proves the theorem. 

Consider now a more complicated problem. Suppose that we have an integral of the form 



X 



dx [f{x) +g{x)]. 



(A12) 



where / and g are both singular on a surface V in the space of the real variables x. Suppose 
that the strength of the singularities are such that the integral of either function would be 
logarithmically divergent. Suppose further that there is a cancellation in the sum such that 
the integral of / + ^f is convergent. Let d{x) be the distance from any point x to the surface 
V. Let us cut out a region of radius R around V and write 



X = lim 



dx f{x) + / dx g{x) 

d>R Jd>R 



(A13) 



Now we wish to explore the consequences of deforming the integration contour for the integral 
of /. Thus we investigate (with the same notation as above) 



I(t) = hm 



dx det I— ^^] f{z{x;t))+ f dx g{x) 

d>R \ dx I Jd>R 



(AM) 



Following the previous proof we find that there is a surface term in the integration by parts 
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— I(t) = lim 
dt R^o 



lim 



tl dx ^iBi^detA^ f 
d<R ^x^' \ " dt \ 



i I dS„ <! Bit detA^ f 



dt 



(A15) 



where the integration is over the surface d = R and dS^ is the surface area differential 
normal to the surface. 

We want to arrange the deformation specified by y'^{x; t) so that dX{t)/dt = 0. For this 
to happen, it is clear that y will have to approach as x approaches the surface V. Then 
Bj^ -^ 5^ and det A ^ 1 as x approaches V. Let the dimensionality of the singular surface V 
be N — a. If the function / was such that the original integral was logarithmically divergent, 
then / oc R'"^ for R —^ 0. The integration over the surface gives a factor dS^ oc R"-'^ for 
i? — s> 0. Suppose that the deformation vanishes proportionally to R''. Then 



-—Tit) oc lim 
dt ^ ' R-^o 



J^a-lj^aj^b 



(A16) 



Then dX{t)/dt = if 6 > 1. The choice made in the main text of the paper is 6 = 2. 



35 



REFERENCES 

[1] D. E. Soper, Pliys. Rev. Lett. 81, 2638 (1998). 

[2] R. K. Ellis, D. A. Ross, and A. E. Terrano, Nucl. Pliys. B178, 421 (1981). 

[3] D. E. Soper, beowulf Version 1.0, |http://zebu.uoregon.edu/~soper/beowulf/| . 

[4] Z. Kunszt and D. E. Soper, Pliys. Rev. D 46, 192 (1992). 

[5] G. Sterman, Pliys. Rev. D 17, 2773, 2789 (1978). 

[6] D. E. Soper, Beowulf 1.0 Technical Notes, [http://zebu.uoregon.edu/ "soper/ beowult"/ 



[7] Z. Kunszt, P. Nason, G. Marcliesini and B. R. Webber in Z Physics at LEPl, Vol. 1, 
edited by B. Altarelli, R. Kleiss ad C. Verzegnassi (CERN, Geneva, 1989), p. 373 

[8] Yu. L. Doksliitzer, contribution to the Workshop on Jets at LEP and HERA, Durham, 
December 1990; S. Bethke, Z. Kunszt, D. E. Soper and W. J. Stirhng, Nucl. Phys. B370, 
310 (1992), Eri&tum-ibzd. B523, 681 (1998). 



36 



